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Summary 


The nonequilibrium flow field resulting from the turbulent mixing 
and combustion of a supersonic axi symmetric hydrogen jet in a supersonic 
parallel ambient coflowing air stream is numerically analyzed. The effective 
turbulent transport properties are determined by means of the (k-e) model of 
turbulence (two-equation model) . The finite-rate chemistry model considers 
eight elementary reactions between six chemical species, H, 0, H 2 O, OH, O 2 
and H2- The governing set of nonlinear partial differential equations is 
solved by an implicit finite- difference procedure. 

Radial distribution at two downstream locations of some important 
variables affecting the flow development, such as the turbulent kinetic 
energy, turbulent dissipation rate, turbulent scale length and viscosity, 
are obtained. The results show that these variables attain their peak values 
at the axis of symmetry. The computed distributions of velocity, temperature 
and mass fractions of the chemical . species present give a complete description 
of the flow field considered. 

A possible direct analytical approach to account for species- 
concentration fluctuations on the mean production rate of these species (the 
phenomenon of unmixedness) is also presented. However, the use of the method 
does not seem justified in view of the excessive computer time required to 
solve the resulting systems of equations. Comparison of species mass fractions 
at a distance downstream equal to eight hydrogen jet diameters without unmixed- 
ness and with it, assuming that local- equilibrium turbulence prevails in the 
flow field, revealed no appreciable differences' in the predicted values of 
species mass fractions at these distances. 
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Notation 


backward reaction rate constants [m J /kmol*s, 

6 2 

m /kmol • s] 

specific heat at constant pressure [j/kg*K] 
production rate of species i in reaction j [kmol/nr • s] 

o 

concentration of species [kmol/iir ] 

production rate of species K due to all chemical reactions 

7 

concentration of catalyst C = c ^ [ktnol/m^] 

i=l 

constant coefficients appearing in turbulence model 
forward reaction rate constants [m /kmol's] 
partial static enthalpy of species i [j/kg] 
total (stagnation) enthalpy [m /s ] 
specific turbulent kinetic energy 
dissipation length scale 
Lewis number 

molecular weight [kg/kmol] 
static pressure 
Prandtl number 

radius measured from axis of jet 

-1 -I 

universal gas constant [J kmol -K ] 
production rate of species i due to all chemical 
reactions divided by the total density, = c^/p, 

Eq. (12) 

Schmidt number 
time 

temperature ( °K) 


u 

V 


X 


X 


Y i 



Subscripts 

C 

D 

E 

i 


I 

K 

t 


x- component of velocity 

r- component of velocity 

distance along central axis 

■unknown variables in the model equation (l4) 

catalyst 

mass fraction of species i divided by its molecular 

weight, Y i = c \/\ [kmol of i/kg^J 

mass fraction of elemental hydrogen divided by the 

molecular weight of atomic hydrogen 

mass fraction of elemental oxygen divided by the 

molecular weight of atomic oxygen 

second- order correlations of species concentration 

fluctuations 

mass fractions of species i 
turbulent dissipation rate 
coefficient of viscosity 
density 

constant coefficients appearing in the turbulence model 


chemistry 

diffusion 

external edge of flow 

value of quantity at initial station; number of species 

i = 1, .... 7 

internal edge of flow 

number of species^ K = 1, . . . , 7. 

refers to turbulent quantity 



r 



nondlmensional quantity 
time- averaged value 


« 
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1 . INTRODUCTION 


Chemical reactions occurring in turbulent flows are of great 
importance in several engineering problems, e.g., in industrial furnaces, 
air-breathing propulsion systems, rocket exhaust plumes, chemical lasers 
and dispersion of pollutants introduced into the atmosphere. The investi- 
gation of such flows is complicated by the simultaneous occurrence of fluid- 
mechanical and chemical effects, the reactants first mixing and then combining 
so as to release chemical energy. A couplet e analysis of such processes is 
a formidable task. Consequently, many simplifications are required for a 
solution and insight into the problem. As in all viscous flows involving 
chemical reactions, it is appropriate to consider two characteristic times, 
one associated with the turbulent diffusion process t& and one with the 
chemical reaction t r - If t^ » t r ? the reaction process is much faster than 
the turbulent diffusion process. Then the flow may be assumed to be in 
chemical equilibrium everywhere. However, in most practical situations ta 
and t r can be of the same order of magnitude and, therefore, the interplay 
of chemical and fluid dynamic effects must be considered. In such instances, 
calculations of the flow are made under the assumption that the mean rate of 
production of each species equals the instantaneous rate with the mean 
quantities such as temperature, density, and species concentration replacing 
their instantaneous values. This assunption does not account for the effects 
of turbulence on the reaction rate which may be significant in certain regions 
of the flow field. 

In recent years the problem of free- jet turbulent mixing coupled 
with chemical reactions has attracted the attention of several groups of 
investigators (Refs. 1-12). The problem considered in all of these investi- 
gations is basically the same, viz. the turbulent mixing and combustion of a 
plane or axisymmetric straight primary fuel jet (or rocket exhaust plume) in 
a parallel ambient coflowing medium. Different investigators have approached 
the problem by introducing different simplifying assurptions about the 
turbulence model describing the turbulent-mixing process, and the chemical 
model describing the particular combustion process considered. Some of the 
authors have assumed that the flow field is in local chemical equilibrium, 
while others have considered the nonequilibrium effects which are determined 
by the chemical kinetics of the reactions involved. However, in almost all 
the investigations the turbulent mixing process was treated by simple 
"algebraic" turbulent-viscosity models or by "differential" turbulent- 
viscosity models. The latter involved one differential equation for the 
turbulent kinetic energy. The other property of turbulent motion required 
to determine the turbulent viscosity, viz. the length scale of turbulence, 
was obtained from various assumed algebraic relations. 

In the present investigation the problem of turbulent mixing and 
combustion of a round, straight hydrogen jet in a parallel ambient coflowing 
air is treated by modelling: (a) the turbulent mixing process with a system 

of two partial differential transport equations for the turbulent kinetic 
energy and its dissipation rate which define the eddy viscosity, and (b) the 
hydrogen-air combustion process with a set of differential equations describing 
the finite- rate production of the species considered. For the set of reactions 
considered, an analytical- numerical approach, based on the use of second-order 
correlations of the species concentrations and their balance equations properly 
closed in accordance with the turbulence model used, is suggested to account 



for the effects of turbulent fluctuations of species concentrations on the 
mean species production rate. An improved version of the CHARNAL code 
(Calculator of Hydrogen-Air Reactions for NASA Langley) (Ref. 13) developed 
at NASA Langley Research Center, was used extensively in the solution of 
the problem. 

Radial distributions are obtained at two downstream locations of 
some important variables affecting the flow development such as the turbulent 
kinetic energy, turbulent dissipation rate, turbulent length scale and 
turbulent viscosity. The conputed distributions of velocity, temperature 
and mass fractions of chemical species give a complete description of the 
flow field considered. The numerical predictions are compared with two sets 
of experimental data. Good qualitative agreement is obtained. 

The continued interest shown by Professor I. I. Glass and his 
constructive comments throughout the course of this work is greatly 
appreciated. Thanks are due Dr. H. L. Beach, Jr., NASA Langley Research 
Center, for providing the opportunity to carry out the present research. 

The valuable and considerable help and cooperation obtained from Dr. J. S. 
Evans, NASA Langley Research Center, in performing the numerical computations 
throughout the course of this investigation is gratefully acknowledged. It 
is a pleasure to thank the Director, Professor J. H. de Leeuw, and the Staff 
of the Institute for Aerospace Studies, University of Toronto, for providing 
me with this opportunity to work at UTIAS. The assistance of Mrs. Winifred 
Dillon in preparing the manuscript is sincerely appreciated. 


2. COMBUSTION OF HYDROGEN IN AIR AT CONSTANT PRESSURE 


The importance of finite- rate chemistry in accurately predicting 
the flow field properties has been shown in various studies (see, for example, 
Refs. 1, 6, 7, 8). The accuracy of the predictions depends largely on the 
employed reaction mechanism. It should be noted, however, that each additional 
species adds to the computational difficulties associated with solving the 
differential equations in the mathematical model of the reacting system. The 
reaction mechanism for hydrogen-air combustion is fairly well established 
(Refs. 14-18). The rate constants for all reactions in the chain have 
been determined but the accuracy of these determinations can be questioned 
to possibly within an order of magnitude. Any inaccuracies in the rate 
constants are naturally carried on to the temperature and composition 
histories which result from a detailed analysis of the system. It has been 
found that eight forward and eight backward reactions involving six reacting 
species plus nitrogen as an inert gas* describe the essential mechanism at 
temperatures of interest. These are listed as follows: 


1. H + 0 = OH + 0 

b i 


2. 0 + H p = OH + H 

b 2 


3. Ho + OH = H + H p 0 

b 3 


4 . 2 OH = 0 + H 0 

b 4 


* For temperatures below about 2200°K the nitrogen reactions and subsequent 
formation of NO and N have been shown to be unimportant (Ref. 17) . 
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5. Hp + X == 2H + X 
b 5 


7 . 0H+X=0+H+X 

b„ 


6 . H^O + X = OH+H+X 


*8 

8 . 0 + X = 2 0 + X 

b 8 


(1) 


where X is a catalyst. The six active species H, 0, H 2 O, OH, O 2 , H 2 and the 
inert species N 2 considered are numbered respectively from 1 to 7* The 
corresponding reaction rate coefficients were estimated as follows: 


3 * lo^e- 8810 / 1 

b l 

= 2.48 x 10 13 e" 660//T 

3 x lO 14 *- 1 * 030 / 1 

b 2 

= 1.30 x 10 l4 e- 24 90 /T 

3 x io 14 .- 3 ^ 0 / 1 

b 3 

= 1.33 x 10 15 e- 10 95 °A 

3 x ^V 3020 / 1 

b 4 

= 3.12 x 10 15 e- 1251 °/ T 

18.5 x 10 ^-i e -^000/T 

b 5 

= 10 l6 

9.66 x lO^VV 62200 / 1 

b 6 

= io 17 

8.0 x ^VV 52200 / 1 

b 7 

= 10 l6 

5.8 x ^- 1 ^ 60600 /T 

b 8 

= 6 x 10 l4 


(2) 


3 

f 4 
f 5 
f 6 
f 7 
f 8 

where T is the absolute temperature in degrees Kelvin. The rates fp-fQ 
and bp-blj. are in units m3/kmol.s, while bj-bS are in mP/lonol2 .s . 

Let ci be the concentration of species i (in kmol/m3) . The 
concentration C of the catalyst X is assumed to be approximately the sum of 
the concentrations of all species including N2: 

i=l 

Then, at any instant of time, the rate of change of the concentration of 

species i in reaction j, c. ., is: 

a J 
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d ll “ ‘ f l c l c 5 + b l c 2 c 4 

d 36 = 

” d l6 

d 12 = f 2 C 2 C 6 ' Vi c 4 

d 4l = 

d 21 = ~ 6 11 

d l3 = f 3°4 C 6 ~ b 3°l C 3 

d 42 = 

C 12 

d l5 = 2 f 5 c 6 C ~ 2 Vl 2c 

d 43 = 

' C 13 

d l6 = f 6 c 3 C - b 6 c 4 c l C 

d 44 

- 2d 24 

d l7 = f 7 C 4 C ‘ Vl C 2 C 

d 46 = 

d l6 

C 21 = ' C ll 

c 47 = 

" c 17 

C 22 = -C 12 

d 51 = 

d ll 

c 24 = f 4 c 4 “ b 4 c 2 c 3 

6 58 = 

"W 2 

c 27 = C 17 

d 62 = 

" C 12 

5 28 ■ 2 f 8 c 5 ° - “a^ 0 

d 63 = 

" C 13 

c 33 = d 13 

d 65 = 

“ d l 5 / 2 


(3a) 


c 34 c 2 4 

The total rate of change of the concentration of species i due to 
all of the chemical reactions is: 


6 1 = d ll + d 12 + d 13 + d 15 + d l6 + d l7 
d 2 = “ d ll ' 6 12 + d 24 + d 17 + d l8 
d 3 = d l3 + d 24 " d l6 

d 4 = " d ll + d l2 ' d l3 " 2d 24 + d l6 ' d l7 


(3) 

Contd. . . 
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( 3 ) 


S = 6 H ' W 2 

6 6 = " (6 12 + 6 l 3 + 6 15 /2) 


The conservation of the number of atoms of hydrogen and oxygen 
(element conservation equations) are obtained respectively as linear com- 
binations of Eq. (3) 


c 1 + 2c^ + c^ + 2cg = 0 

(4) 

c 2 + + 6^ + 26 ^ = 0 

The concentration of species c. are related to their mass fractions 
a ^ as follows: 


c i 



(5) 


where is the molecular weight of species i, and p the total density of 
the gas mixture . 

The conservation equations for a complex reacting gas are (see, 
for example, the book by Penner, Ref. 19) : 


Conservation of species: 


dY i 

P i— 1,...,6 (6) 


Conservation of energy: 



(7) 


Equation of state: 


p = RpYT 


( 8 ) 


Hydrogen element conservation: 


Y f - 


Y. + 2Y_ + Y,. + 2Y, 


(9) 


5 



Oxygen element conservation: 


Y 4 + 2Y 5 (10) 

7 

= Y = ^Y^, and is the 

i=l 

partial specific static enthalpy of species i. We will assume that the 
reaction takes place at constant pressure. Using Eqs. (9) an d (10) to 
determine Y5 and Yg, we can reduce the number of equations in (6) to four 
and, performing the differentiation in Eq. (7), we get the following set 
of six equations to determine the four species Y^, Y^, Y^, Y^, the total 
density p and the temperature T: 


Y x - Y 2 + Y 3 + 


where R is the universal gas constant, Y^ 



i = 1, . . . , 4 


dT 

dt 


I K T ) 3 i 


i=l 




i=l 


( 11 ) 


( 12 ) 


dp = _ £ dl p dY 
dt T dt ' Y dt 


(13) 


•X" 

where Si = Ci/p and hi = hiMi/RT. The nondimensional enthalpy h* is calculated 
by the polynomials 


* -A^i A3i 2 ^4i 3 ^5i k 

h i = A ii + T + T + nr T + T + V 1 4 


if T < 1000° 


and 


* B 2i B 3i 2 B ^i 3 B 5i 4 B 6i 

h ± = B 1± + — T + T + -jp + ~ T + — = if T > 1000° 


The coefficients of these polynomials were fitted using the computer program 
given in Ref. 36. Thermodynamic properties were obtained from JMAF thermo- 
chemical tables (Ref. 20) . The nondimensional specific heats c^. = d(h* T)/dT . 
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FINITE-RATE CHEMISTRY COMPUTATIONS 


The nonlinear coupled set of differential equations (ll)-(l3) are 
of the form 


dx. 

= s i( x i’ x 2 , xg) (i = 1, 6) (14) 


It is known (see Refs. 21, 22) that the numerical integration of 
these equations poses, in certain cases, a severe numerical stability problem. 
Numerical solutions involving the use of predictor- corrector methods or the 
use of the Runge-Kutta technique can be prohibitively expensive because they 
require excessive computing time. Several attempts have been made to over- 
come the difficulty by developing special numerical methods (Refs. 23-27)- 
In the present investigation the method of Lomax and Bailey (Ref. 26) has 
been selected for the solution of Eqs . (ll)-(l3). The system of equations 
(l4) is solved by the implicit modified Euler method: 


n+1 

x. 


n 

i 


S i +1 + S i 


At 


(15) 


where n is a local reference time level. The value of S^ is calculated 
by expanding each Si uniformly with respect to all dependent vaxiables 
x^ = p, X2 = T, X3 = Yq, xi| = Y2> xc = Y3 and xg = Y4, in a Taylor series 
about the time level n and keeping the linear terms only. 



Substituting Eq. (l6) into Eq. (15) we get 



or in matrix form 

( [I] - |r [E] ) ('x n+1 - 1 c 11 ) = f At (18) 

where [I] is the unit matrix, the elements of the matrix E are eij=^Si/9xj | a , 
and the components of column vectors x and S are respectively xp, ...,xg 
and Si, ...,Sg. For any given set of initial data, the linear algebraic 
system of six equations, Eq. (l8) , can be solved by standard numerical 
techniques. The partial derivatives SSi/Sx^ | n in the matrix E were calculated 
nurneri cally . 
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The technique outlined above has been programed, and test runs 
have been made for the static combustion of a stoichiometric hydrogen-air 
or hydrogen-* Vitiated" air mixture for three sets of values of initial 
pressure and tenperature of the mixture. The rate constants used are 
given in Eq. (2). The results are shown in Figs. 1-7* 

Consider the temperature histories depicted on Fig. 7. These 
temperature histories are characterized by two time periods: the initial 

time interval , called the induction time of combustion, when the temperature 
remains almost constant, and a later period, called the energy release time, 
when the temperature increases rapidly from its initial value to the final 
temperature which is the adiabatic flame temperature. It is seen that under 
conditions of high initial tenperature and relatively low pressure the second 
time interval dominates, while under conditions of relatively low initial 
temperature and high pressure the two time periods are of approximately 
equal duration. 

The changes in composition during the intervals of time are shown 
in Figs. 1-6. The behaviour of the atomic hydrogen is of great thermodynamic 
significance because of its high species enthalpy per unit mass. During the 
induction period the mass fraction of H and the other intermediates 0 and OH 
increase rapidly and approach their equilibrium value. During this time 
molecular oxygen and hydrogen are depleted and water is formed. However, the 
chemical energy released by the formation of water is effectively absorbed 
by the atomic hydrogen so that no significant change in temperature occurs. 

The initial values of pressure and temperature do not seem to affect appreciably 
the maximum values of mass fraction, which are of the order of several percent. 
In the next section the present finite-rate chemistry computation will be 
used in the analysis of turbulent free jet hydrogen-air diffusion flames. 


4. GOVERNING EQUATIONS 


The governing equations for the coaxial mixing of two dissimilar 
gases including heat release by chemical reactions in the mixing region were 
derived assuming that the potential core region and the outer flow region 
are inert, inviscid and uniform. The idealized flow configuration under 
study is presented in Fig. 8. A central axisymmetric jet of hydrogen dis- 
charges into a coflowing airstream serving as the oxidant medium. For 
purposes of the present investigation the time-averaged conservation equations 
are used in a boundary layer form with turbulent Lewis number Let = 1. Assuming 
that turbulent mixing predominates to the extent that molecular transport 
becomes negligible, the governing equations in the axisymmetric x, r coordinate 
system take the following form: 

Continuity: 


57 ( pu ) + 7 


d(rpv) = 0 

or 


Axial Momentum: 




dx 



r 


8u 

3 ? 


) 


8 


(19) 

Contd. 



II 


Species Consei r vation: 


Sy 


pu 


K 
cix 


Sy 


1 

r &r 


dY. 


) + <V(P> T > Y i> y 


Sc, 




K 


3’ 


V 


(K = 1, 


4 ) 


Hydrogen Element Conservation: 

SY SY f lh & f X 

p U^ r + pv^ r --^;^^- r ^ J 

Energy: 



SH 1 
pV ^ = r 


d 

& 



(H - k) 


1 



r S(u 2 /a) 

dr 


(19) 


where all the dependent variables involved represent time-averaged flow 
quantities* u and v being the velocity components in the x and r directions, 
respectively, H the total enthalpy of the fluid element, and Pr-fc and Sc-fc the 
turbulent Prandtl and Schmidt numbers, respectively. The quantities Yk 
represent mean species mass fractions divided by their molecular weights: 

Yk = C^k/Mec* In the species conservation equation, Eq. (19) , it is assumed 
that the mean rate of production of each species, cj£, equals the instantaneous 
rate with the mean quantities such as temperature, density, and species con- 
centration replacing their instantaneous values. However, concentration 
measurements in turbulent flames with highly active reactions show that the 
time-average reaction rate in such flames is less than the values given by 
using the time-average values of the species mass fraction. This effect is 
due to high local fluctuations of temperature and species concentrations 
that may give rise to situations where the oxidizer and the fuel are not at 
the same place at the same time and may cause variations in the Arrhenius 
reaction rates, thus reducing the average reaction rate to values lower than 
those which would exist if the local average concentrations and temperature 
would prevail for the entire period. In the literature, this phenomenon is 
called "unmixedness" . The analytical investigation of the combined effect 
of temperature and concentration fluctuations on the mean rate of production 
of species for a realistic combustion process in turbulent flows is a for- 
midable task. Borghi (Ref. 28) has considered this problem for the simple 
case of a single forward bimolecular reaction K + vO -»P, where K denotes 
the fuel, 0 the oxidizer, P the reaction products and v is the stoichiometric 
coefficient, assuming small temperature fluctuations. Later Borghi et al 
(Ref. 29 ) used approximate forms of probability density distributions for 
species production rate. Recently, Spiegler et al (Ref. 30) presented an 
interesting simplified model of unmixedness and applied it to the case of 
hydrogen-air axisymmetric diffusion flame. The results do show some improve- 
ment in the agreement of predicted concentration profiles with experimental 
results. A possible direct analytical approach to account for species 
concentration fluctuations only is presented in Appendix A. 

The following initial and boundary conditions are imposed on the 
above parabolic partial differential equations, Eq. (19): at the initial 

section, x = 0, velocity, temperature (or total enthalpy), species mass 
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fraction profiles should he given. If possible , these initial profile shapes 
should be based on experimental data for the particular configuration of 
interest. In the absence of experimental data, suitable approximations for 
these profiles can usually be provided on the basis of general knowledge of 
two-dimensional or axisymmetric boundary-layer and channel flows. At the 
axis of symmetry the boundary conditions are: 


r 


du _ dH _ _ 


(20) 


and. at the outer edge of the mixing region the dependent variables should 
tend to their values in the outer (external) inviscid flow: 


r oo, u u T 


H -> 




f -f B , y k -y 


‘k e 


(21) 


5. THE TURBULENCE MODEL 


The eddy viscosity pt in Eq. (19) is determined by means of two 
transport equations for the turbulent kinetic energy, k, and its dissipation 
rate, e (two-equation model of turbulence) (Ref. 31) . According to this 
model the magnitude of the eddy viscosity depends only on the local values 
of k, € and the fluid density in the following way: 


p, = C 

p 


P k 2 k 3 / 2 

T 5 e " ~T~ 


(22) 


The quantities k and e satisfy the following differential equations: 


<3k 

pu ^ + 


cik 

pv SF 


a 

SF 




+ n+ 



pe 


(23) 


Be Be 1 a 
pU + pv ciF = F c)F 


The values of the constants Cm, 



Ci, C 2 5 and cr G are (see Ref. 13): 


(24) 


c 

n 

c_ 

c r 

a ' 

p 

1 

2 

k 

e 

.09 

1.43 

1.92 

1.0 

1.3 


The values of the Prandtl/Schmidt numbers are (Ref. 13): 


Pr t = Sc t =0.7 (Le t = l) 


(25) 
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6. RESULTS AND DISCUSSION 


The combustion model developed in Sections 2 and 3 to compute 
finite- rate hydrogen- air reactions, together with the ten equations (19 5 
23 and 24) , were used to predict the free turbulent mixing and combustion 
of a hydrogen jet emerging into a co-axial air stream at a matched exit 
pressure of 1 atm. A schematic of the experimental setup of the flow 
configurations considered is depicted in Fig. 8. The air is preheated by 
partial combustion with hydrogen and the oxygen consumed is replaced. The 
composition in mass fractions of this "vitiated air" is given as 0.24l 02? 

0.281 H2O, 0.478 N2; it also contains small amounts of H, 0, OH and H2. 

For a more detailed description of the experimental apparatus, see Refs. 4 
and 32. The numerical technique adopted is that of Ref. 13 5 and an improved 
version of the CHARNAL code (Ref. 13) developed at RASA Langley Research 
Center was expensively used in the solution of the problem. The chemistry 
and the turbulent diffusion processes are uncoupled by the following procedure, 
first devised by Ferri et al (Ref. 33). 

Let 


y k = Y % + Y Kc 


T 




+ T 


C 


Then the species conservation equation (19), may he split into the following 
two equations: 


and 



SlKc 

ay Kc 



pu 

~sr + 

pv -sr = p 

5 dt 

6 k 


5yk d 1 8 

(* t 


pu - 5 F 

+ pv - 

"5y~ r dr 


dr 

equation 

•» (19): 

, splits to 




8 h c 

Bh c 

ffi c _ n 



pu ^T 

+ pv ^r = 

dt 0 



and 


pu 


S V ^ u r l 


dx 


+ pv 


r 


r 5? ( V k) + V 1 - 


_1_ 

Pl\ 


(26) 

(27) 


( 28 ) 


d(u 2 /2) 


r “Sr 


(29) 


Equations (26) and (28) pertain to a premixed flow of constant velocity and 
pressure with finite-rate chemistry and are similar to Eqs. (6) and (7) in 
Section 2, while Eqs. (27) and (29) refer to purely diffusive flows. The 
numerical solution is then performed as follows: 
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Diffusion Step: 

For a given forward stepsize Axp, Eqs. (27) and (29) are used 
along with other conservation laws in Eqs. (19) and Eqs. (23) and (2k) 
to find the solution of the purely diffusive step. 

Chemistry Step: 


The diffusion stepsize Asp is subdivided into a finite number of 
chemical steps. Ate = Axp/nu, where n is the number of subdivisions, and 
the diffusion solution is used as initial conditions for the chemistry step 
which uses Eqs. (26), (28), or equivalently, Eqs. (6)-(l0) to determine, by 
the method presented in Section 3, the final values of the species, density 
and temperature. The value of n is determined by accuracy requirements for 
the chemistry step solution. 


The quantity Y x = OxAfe? where a x is the mass fraction of elemental 
oxygen, is yet to be determined. It satisfies a transport equation which is 
identical to the hydrogen element conservation equation, Eq. (19)- With 
appropriate initial and boundary values the transport equation for Y x can 
be solved and its value determined. However, using the identity of the 
transport equations for hydrogen and oxygen element conservation and intro- 
ducing the nondimensional variables 


Y f - Y fE 
f " YY ic - Iy E 


* _ Y X “ Yx E 
X " Yx ic - Yx E 


(30) 


where subscript E stands for the value of the variable at the external boundary 
of the mixing zone, and the subscript ic for its initial centerline value, which 
take identical boundary values of SY^/Sr = dY^/c)r = 0 and Y? E = Y x ^ = 0, we can 

obtain the same solution from the two identical transport equations for Y^ and 
Y x if we assume identical initial values (subscript i) for Yf and Y* (same 
differential equations, same initial and boundary values) , i.e. 


Yf i ~ Yf E _ Yx i Y * E 
Yf ic" Yf E Yx ic " Yx E 


(31) 


In the jet mixing problem considered, Yf E = af^/lvL, Yf^ c = l/M^, y x e - a x E /M?> 
and. Y x ^ c = 0. Hence from Eqs. (30) and (31) we have 

Y x (1 - M, Yf . ) 

Y % ■ 

It is easy to verify that the initial distributions of hydrogen and oxygen satisfy 
Eq. (32). Consequently, Y* = Y* and from Eq. (30) we get 

w 1 - Vf) 

* 1 - 

The mass fraction of nitrogen is then obtained from the relation 


a f + a x + a 7 = 1 
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Figures 9-20 show some computed results at two downstream locations, 
x/Dj = 6.7 and 15.5, where Dj = 9-525 mm is the diameter of the hydrogen jet 
(Fig. 8). The initial profiles of the quantities described are also shown 
for convenience. With air diffusing into the inner regions and molecular 
hydrogen diffusing from the inner jet core into the outer regions, a turbulent 
mixing zone with chemical reactions is formed. The extent of the mixing 
and combustion process is indicated by the turbulence quantities profiles, and 
velocity, temperature and mass fraction profiles. 

Distributions of some of the important variables affecting the 
flow development are shown in Figs. 9-12. The initial distribution of the 
turbulent kinetic energy is calculated on the basis of the mixing length 
hypothesis, the mixing length being proportional to the thickness of the 
hydrogen nozzle wall d = 1.57 mm (Fig. 8). The quantity d is also used as 
the constant initial turbulence scale (Fig. 11a). The initial distributions 
of the dissipation rate e (Fig. 10a) and the eddy viscosity pt (Fig* 12a) 
are then calculated from Eq. (22) . It is seen from these figures that far 
downstream (curves c) all the quantities involved tend to reach their 
maximum values at the axis of symmetry. As the shear flow develops downstream 
the mixing zone spreads and the levels of the corresponding variables fall. 

At the outer edge of the mixing zone the curves should tend to the same value 
of the corresponding variable which prevail in the outer flow. However, the 
confuted profiles depict this trend only approximately. 

Velocity and tenperature profiles at the selected downstream 
locations are shown on Figs. 13 and ik . The wake region of the profiles 
(Fig. 13 , curves a and b) arises from the wall boundary layers of the hydrogen 
pipe. Note that the potential core (see Fig. 8) is still present at x/Dj = 
6.7. As the mixing progresses downstream, hydrogen jet and air velocities 
tend to equalize. At the outer edge of the mixing zone the velocity and 
tenperature profiles tend to values prevailing in the external stream. From 
the temperature profiles shown in Fig. l4 the region of combustion is clearly 
evident through the peak in the tenperature level. It can be seen that the 
increase in the maximum temperature values becomes less appreciable as the 
mixing process develops downstream. No temperature increase is observed at 
the centreline with a distance 6-7 diameters from the origin. 

As a result of the combustion process, chemical species are formed. 
The calculated distributions of atomic species H and 0 of the hydroxyl radical, 
OH, and water, H 2 O, are shown in Figs. 15-18. It can readily be seen that the 
peak values reached by these species at some distance downstream remain 
approximately constant or decay, and as the shear flow develops and the 
jet spreads, shift toward the outer region of the flow. The maximum computed 
value of OH at x/Dj = 15-5 is ~ 2 by mass, and is about ten times as large 
as the maximum mass fraction of H and approximately three and a half times 
as large as the peak mass fraction of 0. Within the computed distances, 
traces of H are present at the centreline, whereas the mass fractions of 
0 and OH are still very small at the axis. The mass fraction profile of the 
most important combustion product H 2 O displays peaks in the inner region of 
the flame where the exothermic reactions are most intense. Traces of H 2 O 
are observed at the centreline only after x/Lj ^ 2. All profiles tend to 
their values in the outer region of the flow. 
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Og distributions are primarily determined by the diffusion- controlled 
mixing process, but as chemical reactions become important in the turbulent 
mixing layer, the interaction between the diffusional process and the 
chemical reactions become very complicated, as indicated by the peaks and 
minima in the predicted profiles in Fig. 19- The mass fraction distribution 
of molecular hydrogen is depicted in Fig. 20. 

To account for the effects of species concentration fluctuations 
on the mean production rate of species (the unmixedness) , the analytical 
approach developed in Appendix A was programmed and the subroutine incor- 
porated into the main computer program. Local equilibrium turbulence was 
assumed (production of turbulence and its dissipation rate balance each 
other, i.e. CgL = Cg 2 = 0 in Eq. (A13) ) • In view of the excessively large 
amount of computer time required to solve the system of twenty- four differ- 
ential equations with finite-rate chemistry (70 min on CDC 7600 for a distance 
x/Dj = 5) further computations with this method did not seem justified. Com- 
parison of species mass fractions at x/Dg= 8.2 without unmixedness and with 
it-, revealed no appreciable differences in the predicted values of species at 
these distances. 

Analytical predictions were compared to experimental data obtained 
by Cohen and Guile (Ref. 4) and Evans, Schexnayder and Beach (Ref. 3 5). 
Detailed descriptions of the experimental apparatus, instrumentation and 
experimental procedures are given in these references and in Ref. 32. The 
initial velocity and temperature profiles for the two cases considered are 
shown in Fig. 21. The initial mass fractions of species in the outer 
vitiated air flow are: for the Cohen and Guile case, 0.26 02, 0.15 H 2 O, 

0.59 N 2 , and for the Beach case (Ref. 35), 0.24l 02, 0.28l HpO, 0.478 N2- 
Both cases contained also small amounts of H, 0, OH and H 2 . 

The computed and experimental pitot pressure profiles at axial 
distances 5.2 and 8.9 jet diameters for the Cohen and Guile case and at 
6.56 jet diameter for the Beach case are shown in Figs. 23 and 24. In the 
Beach case the theory predicts the minimum pitot pressure and the general 
trend of the spread of the mixing region relatively well. In both cases 
and at all distances, theoretically predicted spread of the mixing region is 
larger than that given by experimental data. However , it is readily seen 
that the slopes of predicted and experimental profiles are almost identical. 
The discrepancy is primarily due to the differences in the predicted and 
experimental minimum values of the pitot pressure. Species mass fractions 
are compared to data in Figs. 24-26. In both cases and at all distances the 
agreement for the inert species N 2 Is good. Quantitative agreement for hy- 
drogen inproves with downstream distance. The agreement for oxygen is better 
in Beach’s case. In both cases a small amount of oxygen is present in the 
hydrogen-rich region near the axis of symmetry. This is due to the diffusion 
of oxygen into the hydrogen region in the gap between the lip of the hydrogen 
nozzle and the point wheiethe flame is initiated. Figure 26 shows that in 
both cases the predicted water profile bumps are not placed according to 
data. Again, relatively better agreement is achieved in Beach's case. 



7. CONCLUSIONS 


By using a two-equation (k-e) model of turbulence and a finite- 
rate chemistry model , an analytical-numerical investigation was made of the 
turbulent mixing and combustion of a round hydrogen jet in a parallel ambient 
coflowing air stream. Theoretical calculations were unable to correctly 
predict experimental results. An attempt was made to account for the 
species concentration fluctuations on the mean-production rate of these 
species (the phenomenon of unmixedness) . Comparisons of species mass 
fractions without unmixedness and with it, assuming that local equilibrium 
turbulence prevails in the flow field, revealed no appreciable differences 
in the predicted values of species mass fractions. Inproved experimental 
data, more accurate determination of initial conditions, preferably given by 
the solution of the full Navier-Stokes equation in the near flow field, and 
proper account of the phenomenon of unmixedness and compressibility effects, 
could considerably decrease the discrepancy between predicted and experimental 
values . 
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APPENDIX A 


DERIVATION OF THE DIFFERENTIAL EQUATIONS 
FOR THE TRANSPORT OE SECOND- ORDER CORRELATIONS 
OP SPECIES CONCENTRATION FLUCTUATIONS 


1 . Species Production Rate in Turbulent Flow 

Let Yk be split, in conventional fashion, into time-average 
parts Yk and. fluctuating part Y£. Thus 


Then, substituting expressions (5) and (Al) into the rate of change of the 
concentration of species i in reaction j, given by Eq_. (3a), and taking the 
average, we obtain the following expressions (temperature fluctuations are 
neglected) : 

*u = p2 [- f i (? A - z i 5 ) + b i (I A + Wi 

*12 “ P 2 [*2<*2*6 + Z 26 > - VW + 

*13 = p2 [ f 3 ( V6 + V - + Z 13 n 

*15 = P 2 [ 2f 5 <V + %)] - p 3 [ 2 t> 5 (\ 2 i + 5 z n + 2i 1 z 1 i ] 

(A2) 

+ z 3 )] - p3 t b 6 ( Vi f + 2z ik + \ z k + Vi 5 ! 

6 17 = P 2 t f 7 (i ¥ + z 4 )] " e 3 tVW + fZ i2 + Vi + L Z 2^ 

“c 24 = P 2 [f 4 (Y 4 2 + z 44 ) - b 4 (Y 2 ? 3 + z 23 )] 

d 28 = p 2 [2fg(Y Y 5 + Z 5 )] - p 3 [2bg(Y Y 2 2 + Y Z 22 + 2Y 2 Z 2 )] 

where 

Z K“^-i Z lK (a, f>, K = l. •••. 7) (A3) 

i=l 
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and 



from the equation of state, Eq. (8), assuming that 

i=l 

Hence, in Eq. (A2) we have the following set of unknown second-order corre- 
lations of species mole-mass fraction fluctuations: 


Z 11 Z 12 

Z 13 

z i4 ; 

Z 15 

Z l6 

Z 17 


Z 22 

Z 23 

Z 24 ■ 

Z 25 

z 26 

Z 27 



Z 33 

z 3 4 • 

Z 35 

Z 36 

Z 37 

(a4) 



z 44 ; 

Z 45 

z 46 

Z 47 





Z 55 

Z 56 

Z 57 




• 


Z 66 

Z 67 


In the next section it 

quantities 

will 

be shown that the following 15 

independent 

z ff 

z fl 

Z f2 

Z f3 

Z f4 




z ll 

Z 12 

Z 13 

z i4 





Z 22 

Z 23 

z 24 


(A5) 




Z 33 

z 3 4 







Z 44 




can he determined from 15 transport equations derived by using the four species 
conservation equations and the hydrogen element conservation equations for the 
instantaneous flow. The quantities on the right-hand side of the dotted line 
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in Eq. (a 4) are functions of quantities in Eq. (A5) and can be determined 
by using the hydrogen and oxygen element conservation equations, Eqs . (9) 
and (10) , in the following way: From Eqs . ( 9)5 (10), (32), (33) and (Al) 
we get: 


Y x = q>(l - M,Y f ) 

h ■ - * 2 - h - i ‘. )/2 

*6 ' <*f - 

i 7 = (1 - q>)(l - MjYjJ/lty 

for the mean values, and 

Y x = " ^ Y f = CT f 

Y 5 = (Y i - Y 2 ' Y 3 " Y 4 )/2 

y 6 = (Y f ' Y l " 2Y 3 _ Y i )/2 

= - (1 - cp)M L Y^./M 7 = BY^, 


(a6) 


(A7) 


for the fluctuating parts. Multiplying Y£, Yg, Y| successively b y Y^., Y^, 
Y^, Y^ ? Y^, Y£, Yg and Y| and taking the average we get: 

Z f5 = ^ cz ff ” z f2 " z f3 " Z f4^/ 2 
Z f6 = ( Z ff " Z fi " 2Z f3 ” Z f4^/ 2 


Z f7 BZ ff 


Z 15 ( cz fl " Z 12 ” Z 13 ' Z l4^ 2 
^25 — ^ Z f2 ~ Z 22 - Z 23 ” Z 24 ^ ^ 2 


Z 35 ^ CZ f3 Z 23 " Z 33 " Z 34^ 2 


(A8) 

Contd. 
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r 


Z W . (CZ fl4 - Z 24 - Z^ - Z^J/2 
Z 55 - (C% f + Z 22 + Z 33 + \ h - 2CZ f2 - 2CZ f3 - 
- 2CZ f4 + 2Z 23 + 2Z 24 + 2Z 34 )/4 

Z l6 = ^ z fl " z ll ' 2Z 13 " Z l4^/ 2 

Z 26 = ^ Z f2 " Z 12 ” 2Z 23 “ Z 24)/ 2 

z 36 = ^ Z f3" Z l3 " 2Z 33 ” Z 34^/ 2 

z 46 = ^ Z f4 " Z l4 ” 2Z 34 " z 44 ) / 2 

z 5 6 = t CZ ff - Z f2 ' < 2C + Vhs ~ < c + ^ Z f4 - CZ fi + 

+ z !2 + z i3 + Z l4 + 2Z 23 + 2Z 33 + Z 24 + 3Z 34 + z 44^ 

z 66 = ^ Z ff + Z ll + ^ Z 33 + Z 44 ” 2Z fi ” i+z f3 ” 2Z f4 + 

+ 4z 13 + 2z l4 + 4 z 34 )/4 

z 17 = BZ fl’ Z 27 = BZ f2 ; Z 37 = BZ f3’ Z 47 = BZ f4 

Z 5Y — ( CZ f7 - Z 27 ” Z 37 - z 47^/ 2 

Z 67 = (Z f7 " Z 17 - 2Z 37 ‘ Z 47 } / 2 (a8) 


2 . Equation for the Transport of the Second- Order Species Concentration 
Fluctuations 


The species and element conservation laws for the instantaneous flow 
can be expressed by the following single equation: 


a 


( pu V + 1 1? < rpv V 


i a 

r ^r 


r sssr 


+ c, 


a 


kmol of a 


nr 


(A9) 


where the subscript C t denotes any element (in that case = 0) or species 
present. Decomposing the dependent variables in Eq. (A9) into their time- 
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average and fluctuation pants , multiplying by Yg and using the continuity 
equation for the mean and fluctuating parts 3 after some transformations we 
get: 

rYV - rW 1 * 

<»“>' ^ -S + * Y e sr + (<”)’ y b sr + Y <i y b + 


at- _ ar ’ y - n y > v 

+ (pv)' y^ + pvy- ^-+#^[r(pv)' =/|: 


n 1 
r sc "S?” 


i a 

r Sr 


M_ "a 

sc ^ 


bY' n bYl, bY' 


n a . P 

-r— s I" I A CL, t X ' qn 

Sc or or pa pa 


(A10) 


Interchanging a and p and adding the resulting expression to Eq. (10) we 
obtain: 

+ ^ ^6 + (pu) . Y . ^ + (p »). ^ ^ + |- [M'ynt'] - 


~5F 


x p sr 


a dx "5c 


Y* 


5_ 

„ 1 

Y aS 

, ^ 1 

l s ; 

r .. u i 

Sr 

r s5S^ 

+ T s? 

Sc 5r" _ 

+ r Sr 

1 

\ 3 
Jra 


5b 


oY c 


Sy ' 5y ' 


-iitr(pv)' Y^Y’] - (pv)' Y' (PV)’ Y ’ ^ - 2 ^ ^ + 

(All) 


+ Y p 6 a + Y a d p + Y p 6 a + Y a 6 p 


Neglecting the turbulent diffusion terns (pu)' Yg and (pu) T Y& in the x- 
direction (in the boundary layer approximation) and the third- order corre- 
lations 3 and averaging we get: 


- ai i Y e ^ ^ 1 8 


pu 


dx 


+ pv 


S r r "§r 


Sy ! Y 1 


5Y SY B ,, 

- W' y pF' ( pv ) ,Y a^- 2 fe^r^? +Y p 6 a +Y a d p (A12 ) 


Using the approximations 

■ ss; 

and modelling the dissipation term (see Ref. 3*0 by: 


x a x p Sc t “Sr 5 Y P Sc fc Sr -5 


p t 3Y a 
- ( pv) 1 Y 1 = -2- 
VP ; a Sc t ck 7 ^ 
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5y' <3y' 

P tL “ P = c SI y-y. 

^ sc ST °g 2 k arp 

and using Eq. (A3) we finally obtain the transport equations for the quantities 

Z ap : 


- Sz ap . - Sz ap l 5 

mi + pv — c-- 2 - = - 


pu 


"ST 


dr 


5? 


p, \ dZ 

H_ + __L ) ^ 

Sc Sc^. J dr 


'op 


m a h2 B oe 

+ °g 1 ^t ^r" ‘ 'Sr ' C g 2 r Z op + Y p c a + Y a c p 


(A13) 


where Cgp = 2/Sct* For the case considered in the present study, a, p take 
the values of f (for the hydrogen element), 1, 2, 3 5 *+ for H, 0, H 2 O and OH, 

respectively Hence, Eq (A13) constitutes a system of 15 equations for the 

15 unknown quantities represented in Eq. (A5). It is assumed that the 
turbulent Schmidt number Set = 0.7 and the constant Cg 2 = 2-0 (Ref 3*0 in 
all the 15 equations. The last two terms on the right-hand side of Eq. (A13) 
can be calculated by developing the Cj_j given in Eq. (3a) into their mean and 
fluctuating parts: 

5 11 + 6 il ■ o 2 t- f l( J l S 5 + Y i Y 5 + Y 5 Y 1 + Y i Y 5> + b l< S A + Y A + Y i ? 2 + Y 2 Y i» 

^12 + q 2 - P 2 [ f 2< Y £ Y 6 + Y 2 Y 6 + Y &> * Y 2 Y 6 > - b 2< Y A + Y A + Y 4 Y 1 + Y A” 

\ 3 * 4 13 ' p2 t f 3 (Y ». Y 6 + y 4 y 6 + Y A + Y i Y 6> ' b 3 (Y l Y 3 + Y i Y 3 + Y 3 Y 1 + Y i Y 3 ,] 

i 15 + q 5 = p 2 [2f 5 (S 6 ! + up + rs 6 + Y’Y 6 ) ] - P 3 [2b 5 (i 1 2 Y + 25 ^ + yyj 2 + 

+ Y 2 Y' + SY-jY-jY' + Y^Y ' ) ] 

\ 6 + q 6 = P 2 [f 6 (^ + YY’ + Y 3 Y' + Y'Yp] - p 3 [b 6 (Y 1+ Y 1 Y + Y^ + YY 4 Y* + 

+ YY£Y 4 + Y^Y’ + Y-^Y ' Y 4 + Y 4 Y'Y^ + Y-jY^Y’)] 

\ 7 + = p 2 [f ? (YY 4 + YY4 + Y'Y 4 + Y'Yjp] - p 3 [b ? (Y 1 Y 2 Y + YgYY^ + Y^Y' + 

+ YY-jYg) + Y^Y’ + YgY^Y 1 + Y^Y’ + Y-pf^Y')] 

i 24 + c’ 4 = p 2 [f 4 (Y 4 2 + 2Y 4 Y4 + Y4 2 ) - b 4 (Y 2 Y 3 + Y^Y 3 + Y^Y 2 + Y^Y’ ) ] 
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*28 + *28 = P 2 [ £f 8 ( ^5 + ^5 + Y ’% + Y ’ Y pl " P 3 [2b 8 (YY 2 2 + 2Y 2 YY' + YY’ 2 + 

+ Y 2 2 Y’ + 2Y 2 Y^Y' + Y^’)] 

7 1 1 7 , 1 1 

using Eq. (3), multiplying the corresponding ca + cq; by Yp and cp + cp by Ya> 
t aking the average and neglecting third and higher order correlations. The 
resulting expressions will then be functions of p, f]_, f8> bp, b8, 

Yj_ and Z^p, and being cumbersome are not reported here. Their derivation is 
straightforward . 
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TIME (sec) 

FIG. 1 VARIATION OF MASS FRACTION OF ATOMIC HYDROGEN a H in 
STOICHIOMETRIC HYDROGEN-AIR MIXTURES AT CONSTANT 
PRESSURE. 

CURVE I: P = 1.51 atm, T = 1340°K 

CURVE II: P = 1 atm, T = 2000°K 

CURVE III: P = 1 atm, T = 1140°K, Vitiated Air <tc, , n . .. 
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FIG. 2 VARIATION OF MASS FRACTION OF ATOMIC OXYGEN IN 
STOICHIOMETRIC HYDROGEN-AIR MIXTURES AT CONSTANT 
PRESSURE. 

CURVE I: P = 1.51 atm, T = 1340°K 

CURVE II: P = 1 atm, T = 2000° K 

CURVE III: P = 1 atm, T = 1140°K, Vitiated Air a 
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FIG. 3 VARIATION OF MASS FRACTION OF WATER a (1 q IN STOICHIOMETRIC 

HYDROGEN-AIR MIXTURES AT CONSTANT PRESSURE. 

CURVE I: P = 1.51 atm, T = 1340°K 

CURVE II: P = 1 atm, T = 2000° K 

CURVE III: P = 1 atrn, T = 1140°K, Vitiated Air a n . . = 

iIqU i m U 




TIME (sec) 

FIG. 4 VARIATION OF MASS FRACTION OF HYDROXYL RADICAL a QH IN 

STOICHIOMETRIC HYDROGEN-AIR MIXTURES AT CONSTANT PRESSURE. 

CURVE I: P = 1.51 atm, T = 1340°K 

CURVE II: P = 1 atm, T = 2000° K 

CURVE III: P = 1 atm, T = 1140°K, Vitiated AirCt - 

lioU imt 





TIME (sec) 

FIG. 5 VARIATION OF MASS FRACTION OF MOLECULAR OXYGEN % 2 IN STOICHIOMETRIC HYDROGEN-AIR 

MIXTURES AT CONSTANT PRESSURE. 

CURVE I: P = 1.51 atm, T = 1340°l< 

CURVE II: P = 1 atm, T = 2000° K 

CURVE III: P = 1 atm, T = 1140°K, Vitiated Air a H9 n j nit = 0.15 






TIME ( sec ) 


FIG. 6 VARIATION OF MASS FRACTION OF MOLECULAR HYDROGEN «h 2 in STOICHIOMETRIC 

HYDROGEM-AIR MIXTURES AT CONSTANT PRESSURE. 

CURVE I: P = 1.51 atm, T = I340°K 

CURVE II: P = 1 atm, T = 2000° K 

CURVE III: P = 1 atm, T = 1140°K, Vitiated Aira^o init = 0-15 




TIME (sec) 

OICHIOMETRIC HYDROGEN-AIR MIXTURES AT CONSTANT PRESSURE. 

n 1 1 


Air i 


= 0.15 






Ox 















FIG. 13 RADIAL VELOCITY PROFILES 
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FIG. 26 RADIAL MASS-FRACTION PROFILES OF WATER 
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